Impact of dissipative effects on the macroscopic evolution of a Vlasov system 
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Numerical diffusion is introduced by any numerical scheme as soon as small scales fluctuations, 
generated during the dynamical evolution of a collisionless plasma, become comparable to the grid 
size. Here we investigate the role of numerical dissipation in collisionless plasma simulations by 
studying the non linear regime of the two stream instability. We show that the long time evolution 
of the Vlasov - Poisson system can be affected by the used algorithm. 

PACS numbers: 52.65.Ff, 52.35.Qz, 52.35.Fp, 52.35.Mw 



Many space and laboratory plasmas can be considered 
as weakly collisional since the collisional frequency is 
smaller than all the other frequencies, as for example 
the plasma frequency. In other words, for these plasmas, 
the mean free path of the particles is (much) longer than 
all the other characteristic length scales of the plasma 
and, sometimes, even larger than the dimension of the 
plasma itself. At first approximations such plasmas can 
be considered as collisionless and their dynamics can be 
well represented using a Hamiltonian description. This 
approach is based on the idea that the dissipative scale 
(for example in numerical simulations the grid size) is 
much smaller than any macroscopic physical length scale 
of the system, so that dissipation has no feedback on the 
macroscopic asymptotic evolution of the system. 

Numerical simulations based on non collisional mod- 
els must necessarily face with the small scales genera- 
tion problem during the dynamical evolution of the sys- 
tem; indeed, when the typical length scales of the fluc- 
tuations become comparable to the grid size, numerical 
dissipation comes into play leading the system to vio- 
late the conservation constraints of Hamiltonian dynam- 
ics and to reconnect close isolines of the distribution func- 
tion (d.f.). This process, formally forbidden, is very well 
highlighted by the time evolutions of the system invari- 
ants Ni = J Pdxdv i = 1,2,.. and by the "entropy" 
S = — J f ln(f)dxdv (here / is the d.f.), showing sudden 
variations when closed vortices are formed in phase space 
as a consequence of particle trapping. 

In this paper, through numerical studies of a non-linear 
regime of a collisionless plasma, we discuss the role of ar- 
tificial dissipation introduced by a numerical scheme on 
the plasma dynamics, influencing the final Vlasov evolu- 
tion of the system even if the grid size is much shorter 
than any physical relevant scale length. The dynami- 
cal non linear evolution we chose for our numerical sim- 
ulations is the well-known two stream instability. In 
this case, dissipation allow for the formation of coher- 
ent macroscopic structures in phase space (vortices). 

Since the non-linear dynamics of the two stream in- 
stability is substantially driven by kinetic effects, es- 
pecially concerning the saturation phase where particle 
trapping play a crucial role, a kinetic approach is nec- 



essary. This can be done using Vlasov equation, which 
replaces Coulomb interactions between charged particles 
with a mean electromagnetic field. This field is deter- 
mined self-consistently trough the particle distribution 
function by Maxwell and Poisson equations. Since the 
two stream instability is driven by purely electrostatic 
mechanisms, we limit our study to the solution of the 
1D-1V Vlasov - Poisson system of equations: 
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In these equations and in the rest of this paper, time 
t is normalized to the inverse of the electron plasma fre- 
quency uj pe , velocities v e and v p to the electronic thermal 
velocity Vth,e, electron and proton distribution functions, 
f e and / p , to the equilibrium particle density no, lengths 
to the Debye length Xo — Vth,e/^pe and the electric field 
E to m e Vth,eU pe I 'e. As mentioned above, the dynamics 
described by this system of equations is characterized by 
the absence of collisions; hence, from a numerical point 
of view, the choice of an algorithm that is capable to con- 
serve better than possible the invariants of the system is 
crucial. Our numerical scheme is based on the splitting 
scheme of Cheng and Knorr, 1976, Q for the solution 
of the Vlasov equation; therefore, the problem is mainly 
reduced to an interpolation problem for the distribution 
function. Here we compare three well known interpola- 
tion algorithms, namely the Van Leer method at second 
and third order Q (at which, in the text, we'll refer as 
VL2 and VL3) and the Spline method Q (a third order 
method). The Poisson equation is solved, at every time 
step, by spectral methods (i.e. fast Fourier Transform 
technique); in particular, we calculate the plasma den- 
sity by integrating the electron distribution functions in 
velocity. 

We made the two stream instability runs for the three 
simulation algorithms (VL2, VL3 and Spline, runs A, B 
and C, respectively) using the same identical parameters 
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FIG. 1: Electron distribution functions in phase space at time t=1200 of, top to bottom, the VL2, VL3 and Spline algorithm, 
runs A, B and C, respectively. 



(however, we recall that the CPU time for the three algo- 
rithms is not equal). The simulations box is L x = 20Ad 
in space and the velocity interval is — v max < v < v max , 
with v m ax = 5v t h,e- We use N x = 128 points in space 



and N v = 501 in velocity, corresponding to a phase space 
grid resolution of dx = 0.16 and dv = 0.02. Other param- 
eters are: amplitude of the initial random perturbation 
of e = 0.0001; dt = 0.0003 w" 1 , total simulation time 
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FIG. 2: The time evolution of the third invariant. The con- 
tinuous, dashed and dot-dashed lines correspond to run A to 
C, respectively. 



equal to 1200 w^ 1 ; modulo of mean velocity of the two 
initial electron streams (uq) equal to 2 v t h.e- The initial 
electron distribution function we use is 



FIG. 3: The time evolution of 8E (same line style as in fig. 
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In Fig. ^ we show the phase space vortices (same con- 
tour levels) generated by the evolution of the instability. 
There is an initial good agreement among the results of 
the three algorithms: the vortex appears at the same time 
for all methods and is displaced in the same position. On 
the other hand, we observe that the Spline vortex (run C) 
propagates with a different velocity with respect to the 
vortex obtained with VL2 and VL3 (runs A, B) method. 
Furthermore, the spatial structure of the VL2 vortex is 
significantly different from that of the vortices obtained 
with the other two algorithms. 

By looking at the behavior of the N% invariant, Fig. 
13 we see in all cases a sudden decrease of the invariant 
as soon as the vortices start to form. This is a conse- 
quence of the d.f. lines reconnection processes at the grid 
scale length where the algorithm is forced to introduce 
some artificial dissipation eventually leading to the clo- 
sure of the particle orbits (i.e. vortices) contrary to the 
Hamiltonian character of the equations. We note that 
both VL3 and Spline invariants tend to a (equal) con- 
stant asymptotic value, while VL2 continues to smoothly 
decrease, meaning that the numerical resolution for the 
VL2 algorithm must be increased (even if both dx and 
dv are much shorter than the vortex dimension), as also 
shown by the time evolution of the total energy varia- 
tions of the system in Fig. [3] where we plot the evo- 
lution of the normalized energy variations defined as 
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FIG. 4: The time evolution of the entropy (same line style as 
in fig. H 



5E = (E tot - E tot {t = 0))/E tot {t = 0). Indeed, we see 
that after the phase space vortex formation, the energy 
variation for VL3 and Spline algorithms becomes nearly 
constant, while for VL2 begins to monotonically increase. 
In Fig. 0]we show the time evolution of the entropy. We 
again observe a strong variation of entropy during the 
vortex formation phase for all three algorithms while in 
the asymptotic limit the entropy becomes nearly constant 
for both VL3 and Spline, while continues to increase for 
VL2. The uncorrect behavior of the VL2 method with 
respect to VL3 and Spline is the consequence of the fact 
that VL2 method is a lower order scheme. We underline 
that VL2 is however a II order scheme (II order schemes 
are often used in collisionlcss simulations) and that the 
grid spacing, dx <C Xd and dv <C Vth,e, seems to be ade- 
quate for correctly describing the formation of a coherent 
structure much larger than dx x dv. To clarify this point, 
we made a number of other runs with different numer- 
ical accuracy, corresponding to different computational 
times. For the algorithms here used, we know that for 
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FIG. 5: The time evolution of the invariants N\ (total charge 
density), N2 and N3. The continuous, dashed and dot-dashed 
lines correspond to runs D, E and F made with the VL2, VL3 
and Spline algorithm, respectively. 



the same grid spacing the computational CPU time scales 
as: 



used is to take a numerical accuracy for the three algo- 
rithms corresponding to the same computational time. 
We made the new runs by using all the same parame- 
ters of the previous runs, but with L x — 30 and v max = 
15^tft,, e - The numerical mesh are: ^=300 and _/V„=601 
for VL2 (run D), ^=200 and JV„=401 for VL3 (run E) 
and N x = 100 and N v = 201 (run F) for Spline. We 
now found that the energy variation is now asymptoti- 
cally constant for all methods. However, even if VL2 has 
a better resolution with respect to the other two meth- 
ods, invariants decrease more than in VL3 or Spline. So, 
even at parity of computational time conditions, the fact 
that VL3 and Spline are third order methods while VL2 
is a second order one has a relevant effect on the invari- 
ants. Finally, VL3 and Spline have different trends for 
invariants, but the same asymptotic values. 

In conclusion, choosing a numerical algorithm means 
to select a determined quantity of artificial dissipation. 
This means that, even if the grid length scales are by far 
the shorter length scales of the system, the final state 
of the system can be affected by the kind of algorithm 
we have used. Furthermore, even by performing very ac- 
curate simulations, i.e. grid scales sufficiently short to 
have a "correct" long time behavior of the energy and 
the invariants, the long time nonlinear dynamics can be 
significantly different, in particular when more "turbu- 
lent" systems are studied (for example if some external 
forcing continues to inject energy during the nonlinear 
regime) . 
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